Hoja de trabajo 03 - Regresión Lineal

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(knitr)
library(ggplot2)
library(cluster)
library(factoextra)

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

require(caret)

## Loading required package: caret

## Loading required package: lattice

library(corrplot)

## corrplot 0.92 loaded

library(Metrics)

## 
## Attaching package: 'Metrics'

## The following objects are masked from 'package:caret':
## 
##     precision, recall

library(boot)

## 
## Attaching package: 'boot'

## The following object is masked from 'package:lattice':
## 
##     melanoma

library(hopkins)
library(e1071)
library(mclust)

## Package 'mclust' version 6.0.0
## Type 'citation("mclust")' for citing this R package in publications.

library(fpc) 
library(NbClust) 
library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

library(FeatureImpCluster)

## Loading required package: data.table

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

library(pheatmap)
library(ggplot2)

options(dplyr.summarise.inform = FALSE)

data <- read.csv("train.csv")

2. Análisis Exploratorio

¿Cuál es la Zona más cara?

grouped_data <- data %>% group_by(MSZoning)
agg_tbl <- grouped_data %>% summarise(median(SalePrice))
grouped_dataset <- as.data.frame(agg_tbl)
colnames(grouped_dataset)[colnames(grouped_dataset) == "MSZoning"] ="Zona"
colnames(grouped_dataset)[colnames(grouped_dataset) == "median(SalePrice)"] ="Promedio de Precio (USD)"
grouped_dataset <- grouped_dataset[order(grouped_dataset$`Promedio de Precio (USD)`, decreasing = TRUE), ]
kable(grouped_dataset, caption = "Promedio de ventas por cada zona")

Promedio de ventas por cada zona
	Zona	Promedio de Precio (USD)
2	FV	205950
4	RL	174000
3	RH	136500
5	RM	120500
1	C (all)	74700

Un 50% de las casas en la zona Floating Village Residential se venden por arriba de los 205940.00 USD. Posteriormente sigue Residential Low Density.

¿En qué vecindario se ubican las casas con mayor área en su terreno?

grouped_area_data <- data %>% group_by(Neighborhood)
area_with_median <- grouped_area_data %>% summarise(median(LotArea))
grouped_dataset_lot <- as.data.frame(area_with_median)
colnames(grouped_dataset_lot)[colnames(grouped_dataset_lot) == "Neighborhood"] ="Vecindario"
colnames(grouped_dataset_lot)[colnames(grouped_dataset_lot) == "median(LotArea)"] ="Area del terreno ft^2"
grouped_dataset_lot <- grouped_dataset_lot[order(grouped_dataset_lot$`Area del terreno ft^2`, decreasing = TRUE), ]
kable(grouped_dataset_lot[0:5,], caption = "Mediana de area por cada vecindario")

Mediana de area por cada vecindario
	Vecindario	Area del terreno ft^2
5	ClearCr	17575.0
24	Timber	12781.5
16	NoRidge	12090.0
17	NridgHt	11694.0
7	Crawfor	11500.0

ggplot(grouped_dataset_lot[0:5,], aes(x=Vecindario, y=`Area del terreno ft^2`)) +
geom_bar(stat="identity", fill="steelblue")

La información nos demuestra que Clear Creek es el vecindario que tiene los terrenos con mayor área. Un 50% de los terrenos tiene un área mayor a 17575.0 pies cuadrados.

¿En qué año hubo más remodelaciones?

data_with_remodelation <- data[data$YearBuilt != data$YearRemodAdd,]

En el dataset hay 696 registros de casas que fueron remodeladas. A continuación agrupamos por año la cantidad de remodelaciones.

grouped_remodelation_data <- data_with_remodelation %>% group_by(YearRemodAdd)
year_with_count <- grouped_remodelation_data %>% summarise(total_count=n(),
            .groups = 'drop')

grouped_remodelation <- as.data.frame(year_with_count)
colnames(grouped_remodelation)[colnames(grouped_remodelation) == "YearRemoAdd"] ="Año"
colnames(grouped_remodelation)[colnames(grouped_remodelation) == "total_count"] ="Cantidad de casas remodeladas"
grouped_remodelation <- grouped_remodelation[order(grouped_remodelation$`Cantidad de casas remodeladas`, decreasing = TRUE), ]
kable(grouped_remodelation[0:5,], caption = "Numero de casas Remodeladas por año")

Numero de casas Remodeladas por año
	YearRemodAdd	Cantidad de casas remodeladas
1	1950	166
44	2006	51
45	2007	42
43	2005	40
38	2000	33

Los 5 años con más casas remodeladas son 1950, 2006, 2007, 2005 y 2000 respectivamente. Es importante notar que en 1950 hay mucha más diferencia que con los demás años.

most_remodelated_zone_1950 <- tail(names(sort(table(data_with_remodelation[data_with_remodelation$YearRemodAdd == "1950",]$MSZoning))), 1)

most_remodelated_zone_2006 <- tail(names(sort(table(data_with_remodelation[data_with_remodelation$YearRemodAdd == "2006",]$MSZoning))), 1)

most_remodelated_zone_2007 <- tail(names(sort(table(data_with_remodelation[data_with_remodelation$YearRemodAdd == "2007",]$MSZoning))), 1)

La mayoría de las remodelaciones pertenecen a la zona RL - Residencial de baja densidad.

¿Cuáles son los materiales más utilizados para la fundación de las casas?

grouped_foundation_data <- data %>% group_by(Foundation)
foundation_with_count <- grouped_foundation_data %>% summarise(total_count=n(),
            .groups = 'drop')

grouped_foundation <- as.data.frame(foundation_with_count)
colnames(grouped_foundation)[colnames(grouped_foundation) == "Foundation"] ="Material para la base"
colnames(grouped_foundation)[colnames(grouped_foundation) == "total_count"] ="Cantidad de casas"
grouped_foundation <- grouped_foundation[order(grouped_foundation$`Cantidad de casas`, decreasing = TRUE), ]
kable(grouped_foundation[0:5,], caption = "Numero de casas por cada material de base")

Numero de casas por cada material de base
	Material para la base	Cantidad de casas
3	PConc	647
2	CBlock	634
1	BrkTil	146
4	Slab	24
5	Stone	6

La mayoría de casas están hechas de concreto. Puede ser que un tipo menos común, pero más caro, como el ladrillo signifique un precio más alto.

¿Qué tipos de casa son más comunes?

grouped_type_data <- data %>% group_by(BldgType)
type_count <- grouped_type_data %>% summarise(total_count=n(), .groups = 'drop')

grouped_type <- as.data.frame(type_count)
colnames(grouped_type)[colnames(grouped_type) == "BldgType"] ="Tipo de casa"
colnames(grouped_type)[colnames(grouped_type) == "total_count"] ="Cantidad"
grouped_type <- grouped_type[order(grouped_type$`Cantidad`, decreasing = TRUE), ]
kable(grouped_type, caption = "Conteo de tipo de casas")

Conteo de tipo de casas
	Tipo de casa	Cantidad
1	1Fam	1220
5	TwnhsE	114
3	Duplex	52
4	Twnhs	43
2	2fmCon	31

Los tipo de casas más comunes son de 1 familia y de tipo town house.

¿Que vecindarios tienen un mayor rango de precios?

grouped_data <- data %>% group_by(Neighborhood) %>% summarise(minPrice = min(SalePrice), maxPrice = max(SalePrice))
grouped <- grouped_data %>% mutate(priceRange = maxPrice - minPrice) %>% arrange(desc(priceRange))

colnames(grouped)[colnames(grouped) == "Neighborhood"] = "Vecindario"
colnames(grouped)[colnames(grouped) == "minPrice"] = "Precio Mínimo (USD)"
colnames(grouped)[colnames(grouped) == "maxPrice"] = "Precio Máximo (USD)"
colnames(grouped)[colnames(grouped) == "priceRange"] = "Rango de Precios (USD)"

kable(head(grouped, 3), caption = "Top 3 - Vecindarios con mayor rango de precios")

Top 3 - Vecindarios con mayor rango de precios
Vecindario	Precio Mínimo (USD)	Precio Máximo (USD)	Rango de Precios (USD)
NoRidge	190000	755000	565000
NridgHt	154000	611657	457657
OldTown	37900	475000	437100

Los vecindarios con mayor rango de precios son NoRidge, NridgHt y OldTown. Se puede observar que el rango de precios de NridgHt y OldTown es similar, sin embargo es necesario tomar en cuenta que el precio mínimo de OldTown es mucho menor al de NridgHt.

¿Cuál proporción de casas en cada zona?

zone_counts <- data %>% count(MSZoning) %>% mutate(Proporcion = prop.table(n) * 100) %>% 
rename(Zona = MSZoning, Cantidad = n) %>% mutate(Proporcion = paste0(round(Proporcion, 2), "%")) %>%
arrange(desc(Cantidad))
kable(zone_counts, caption = "Proporción de casas por zona")

Proporción de casas por zona
Zona	Cantidad	Proporcion
RL	1151	78.84%
RM	218	14.93%
FV	65	4.45%
RH	16	1.1%
C (all)	10	0.68%

Se puede observar que la zona en donde hay mayor proporción de casas es la zona RL, seguido de RM tomando en cuenta que la proporción de casas en esta zona es considerablemente menor así como en las zonas FV, RH y C(all).

¿Cuáles son las características más comunes de las casas que se venden por encima del precio medio?

price_mean_up <- mean(data$SalePrice)
houses_mean_above <- subset(data, SalePrice > price_mean_up)
ggplot(data = houses_mean_above, aes(x = MSZoning)) +geom_bar()

ggplot(data = houses_mean_above, aes(x = BldgType)) +geom_bar()

ggplot(data = houses_mean_above, aes(x = HouseStyle)) +geom_bar()

ggplot(data = houses_mean_above, aes(x = GarageType)) +geom_bar()

Las características más comunes son las siguientes: 1) Zona donde se ubica las casas: RL (Casa residencial de baja densidad) 2) Tipo de vivienda: 1Fam (Unifamiliar) 3) Estilo de la casa: 1Story (Un nivel) 4) Tipo de Garaje: Attchd (Adjuntada a la casa)

¿Cuántas casas tienen piscina o garaje? ¿Como se relacionan estas características con el precio de venta?

print("Counting Pools And Garages")

## [1] "Counting Pools And Garages"

nrow(data) - table(data$PoolArea > 0)

## 
## FALSE  TRUE 
##     7  1453

nrow(data) - nrow(data[data$GarageType=="NA",])

## [1] 1379

print("")

## [1] ""

print("Relation with the seal price")

## [1] "Relation with the seal price"

tapply(data$SalePrice, data$PoolQC > 0, mean)

##     TRUE 
## 288138.6

tapply(data$SalePrice, data$GarageType=="NA", mean)

##    FALSE 
## 185479.5

1453 casas tienen psicina 1379 casas cuentan con garaje

Como se puede ver si influyen la relación con el precio, las piscinas suelen elevar el precio de las casas

¿Cómo se distribuyen las ventas de casas a lo largo del tiempo entre los diferentes periodos?

ggplot(data = data, aes(x = factor(YrSold))) +
  geom_bar(fill = "green") +
  ggtitle("Ventas de casas por año") +
  xlab("Año") +
  ylab("Cantidad de ventas")

ggplot(data = data, aes(x = YrSold, y = SalePrice)) +
  geom_line(color = "green") +
  ggtitle("Tendencia del precio de venta") +
  xlab("Año") +
  ylab("Precio de venta promedio")

Como se puede ver en las dos gráficas presentadas desde el año 2006 al 2009 existio una tendencia donde subía y bajan las ventas de las casas. Pero desde 2010 este fue en bajando.

Preprocesamiento de datos

columns_used <- c()
neighborhoodNames <- c("NoRidge", "NridgHt", "StoneBr", "Timber", "Veenker", "Somerst", "ClearCr", "Crawfor", "CollgCr", "Blmngtn", "Gilbert", "NWAmes", "SawyerW", "Mitchel", "NAmes", "NPkVill", "SWISU", "Blueste", "Sawyer", "OldTown", "Edwards", "BrkSide", "BrDale", "IDOTRR", "MeadowV")

for(n in 1:length(neighborhoodNames)) {
  # Variable minuscula para nuestro uso.
  data$neighborhood[data$Neighborhood == neighborhoodNames[n]] <- n
}
columns_used <- append(columns_used, "neighborhood")

hs <- c("1Story", "2Story", "1.5Fin",   "SLvl", "SFoyer")

for(n in 1:length(hs)) {
  # Variable minuscula para nuestro uso.
  data$houseStyle[data$HouseStyle == hs[n]] <- n
}
columns_used <- append(columns_used, "houseStyle")

 data$houseZone[data$MSZoning == "A"] <- 1
 data$houseZone[data$MSZoning == "C"] <- 2
 data$houseZone[data$MSZoning == "FV"] <- 3
 data$houseZone[data$MSZoning == "I"] <- 4
 data$houseZone[data$MSZoning == "RH"] <- 5
 data$houseZone[data$MSZoning == "RL"] <- 6
 data$houseZone[data$MSZoning == "RP"] <- 7
 data$houseZone[data$MSZoning == "RM"] <- 8
 columns_used <- append(columns_used, "houseZone")

data$houseUtilities[data$Utilities == "AllPub"] <- 1
data$houseUtilities[data$Utilities == "NoSewr"] <- 2
data$houseUtilities[data$Utilities == "NoSeWa"] <- 3
data$houseUtilities[data$Utilities == "ELO"] <- 4
columns_used <- append(columns_used, "houseUtilities")

data$roadAccess[data$Condition1 == "Artery"] <- 1
data$roadAccess[data$Condition1 == "Feedr"] <- 2
data$roadAccess[data$Condition1 == "Norm"] <- 3
data$roadAccess[data$Condition1 == "RRNn"] <- 4
data$roadAccess[data$Condition1 == "RRAn"] <- 5
data$roadAccess[data$Condition1 == "PosN"] <- 6
data$roadAccess[data$Condition1 == "PosA"] <- 7
data$roadAccess[data$Condition1 == "RRNe"] <- 8
data$roadAccess[data$Condition1 == "RRAe"] <- 9
columns_used <- append(columns_used, "roadAccess")

data$remodelated[data$YearBuilt != data$YearRemodAdd] <- 1
data$remodelated[data$YearBuilt == data$YearRemodAdd] <- 0
columns_used <- append(columns_used, "remodelated")

data$roofStyle[data$RoofStyle == "Flat"]  <- 1
data$roofStyle[data$RoofStyle == "Gable"]  <- 2
data$roofStyle[data$RoofStyle == "Gambrel"]  <- 3
data$roofStyle[data$RoofStyle == "Hip"]  <- 4
data$roofStyle[data$RoofStyle == "Mansard"]  <- 5
data$roofStyle[data$RoofStyle == "Shed"]  <- 6
columns_used <- append(columns_used, "roofStyle")

data$roofMaterial[data$RoofMatl == "ClyTile"] <- 1
data$roofMaterial[data$RoofMatl == "CompShg"] <- 2
data$roofMaterial[data$RoofMatl == "Membran"] <- 3
data$roofMaterial[data$RoofMatl == "Metal"] <- 4
data$roofMaterial[data$RoofMatl == "Roll"] <- 5
data$roofMaterial[data$RoofMatl == "Tar&Grv"] <- 6
data$roofMaterial[data$RoofMatl == "WdShake"] <- 7
data$roofMaterial[data$RoofMatl == "WdShngl"] <- 8
columns_used <- append(columns_used, "roofMaterial")

data$overallQuality <- data$OverallQual
columns_used <- append(columns_used, "overallQuality")

data$overallCondition <- data$OverallCond
columns_used <- append(columns_used, "overallCondition")


data$exteriorCondition[data$ExterCond == "Po"] <- 1
data$exteriorCondition[data$ExterCond == "Fa"] <- 2
data$exteriorCondition[data$ExterCond == "TA"] <- 3
data$exteriorCondition[data$ExterCond == "Gd"] <- 4
data$exteriorCondition[data$ExterCond == "Ex"] <- 5
columns_used <- append(columns_used, "exteriorCondition")

data$foundationMaterial[data$Foundation == "BrkTil"] <- 1
data$foundationMaterial[data$Foundation == "CBlock"] <- 2
data$foundationMaterial[data$Foundation == "PConc"] <- 3
data$foundationMaterial[data$Foundation == "Slab"] <- 4
data$foundationMaterial[data$Foundation == "Stone"] <- 5
data$foundationMaterial[data$Foundation == "Wood"] <- 6
columns_used <- append(columns_used, "foundationMaterial")

data$basement[is.na(data$BsmtQual)] <- 0
data$basement[!is.na(data$BsmtQual)] <- 1
columns_used <- append(columns_used, "basement")

data$basementCondition[data$BsmtCond == "Ex"] <- 3
data$basementCondition[data$BsmtCond == "Gd"] <- 2
data$basementCondition[data$BsmtCond != "Ex"] <- 1
data$basementCondition[data$BsmtCond != "Gd"] <- 1
data$basementCondition[is.na(data$BsmtCond)] <- 0
columns_used <- append(columns_used, "basementCondition")

data$fireplace[is.na(data$FireplaceQu)] <- 0
data$fireplace[!is.na(data$FireplaceQu)] <- 1
columns_used <- append(columns_used, "fireplace")

data$garageArea <- data$GarageArea
columns_used <- append(columns_used, "garageArea")

data$pool[is.na(data$PoolQC)] <- 0
data$pool[!is.na(data$PoolQC)] <- 1
columns_used <- append(columns_used, "pool")

data$additionalFeature[is.na(data$MiscFeature)] <- 0
data$additionalFeature[!is.na(data$MiscFeature)] <- 1
columns_used <- append(columns_used, "additionalFeature")

data$livingArea <- data$GrLivArea
columns_used <- append(columns_used, "livingArea")

data$yearBuilt <- data$YearBuilt
columns_used <- append(columns_used, "yearBuilt")


data$salePrice <- data$SalePrice
columns_used <- append(columns_used, "salePrice")

tv <- c("WD", "Oth", "New", "ConLw", "ConLI", "ConLD", "Con", "CWD", "COD")

for(n in 1:length(tv)) {
  # Variable minuscula para nuestro uso.
  data$saleType[data$SaleType == tv[n]] <- n
}
columns_used <- append(columns_used, "saleType")

msz <- c("FV", "RL", "RH", "RM" , "C (all)")

for(n in 1:length(msz)) {
  # Variable minuscula para nuestro uso.
  data$mSZoning[data$MSZoning == msz[n]] <- n
}
columns_used <- append(columns_used, "mSZoning")

Borrando valores inneceesarios

cleanData <- subset(data, select = columns_used)

3. Clustering

Agrupamiento de las columnas más importantes

# Definir las columnas más importantes
important_cols <- c("YrSold", "SalePrice", "YearBuilt", "YearRemodAdd", "PoolArea", "GarageArea", "OverallQual", "OverallCond", "LotFrontage", "LotArea", "GarageCars")
# Normalizar variables numericas
cols_num_norm <- data[,important_cols] <- mutate_if(data[,important_cols], is.numeric, scale)

Resumen de las columnas a utilizar para llevar a cabo el clustering

print(summary(data[,important_cols]))

##       YrSold.V1          SalePrice.V1        YearBuilt.V1    
##  Min.   :-1.3671863   Min.   :-1.838074   Min.   :-3.286697  
##  1st Qu.:-0.6142282   1st Qu.:-0.641296   1st Qu.:-0.571727  
##  Median : 0.1387300   Median :-0.225587   Median : 0.057352  
##  Mean   : 0.0000000   Mean   : 0.000000   Mean   : 0.000000  
##  3rd Qu.: 0.8916881   3rd Qu.: 0.416387   3rd Qu.: 0.951306  
##  Max.   : 1.6446462   Max.   : 7.226343   Max.   : 1.282400  
##                                                              
##    YearRemodAdd.V1        PoolArea.V1        GarageArea.V1   
##  Min.   :-1.6887898   Min.   :-0.068668   Min.   :-2.212205  
##  1st Qu.:-0.8653621   1st Qu.:-0.068668   1st Qu.:-0.647694  
##  Median : 0.4424348   Median :-0.068668   Median : 0.032833  
##  Mean   : 0.0000000   Mean   : 0.000000   Mean   : 0.000000  
##  3rd Qu.: 0.9268040   3rd Qu.:-0.068668   3rd Qu.: 0.481841  
##  Max.   : 1.2174256   Max.   :18.299910   Max.   : 4.420012  
##                                                              
##    OverallQual.V1      OverallCond.V1      LotFrontage.V1       LotArea.V1     
##  Min.   :-3.687150   Min.   :-4.111561   Min.   :-2.01978   Min.   :-0.923413  
##  1st Qu.:-0.794879   1st Qu.:-0.517023   1st Qu.:-0.45502   1st Qu.:-0.296889  
##  Median :-0.071812   Median :-0.517023   Median :-0.04324   Median :-0.104028  
##  Mean   : 0.000000   Mean   : 0.000000   Mean   : 0.00000   Mean   : 0.000000  
##  3rd Qu.: 0.651256   3rd Qu.: 0.381612   3rd Qu.: 0.40972   3rd Qu.: 0.108671  
##  Max.   : 2.820459   Max.   : 3.077516   Max.   :10.00422   Max.   :20.511245  
##                                          NA's   :259                           
##     GarageCars.V1    
##  Min.   :-2.3646297  
##  1st Qu.:-1.0265059  
##  Median : 0.3116179  
##  Mean   : 0.0000000  
##  3rd Qu.: 0.3116179  
##  Max.   : 2.9878655  
##

Hopkins

set.seed(123)
hopkins(data[, important_cols], m=1400)

## Warning in runif(m, min = colmin[i], max = colmax[i]): NAs produced

## [1] 0.9999868

El resultado es de 0.9999868, el cual indica una tendencia a clustering alta ya que es un valor cerca a 1.

Evaluación Visual de Tendencia de la data (VAT)

data_dist <- dist(data[, important_cols])
fviz_dist(data_dist, show_labels=F)

Como se puede observar el método de hopkins esta en lo correcto y se puede confiar en el resultado dado de VAT ya que muestra patrones de agrupamiento en general.

Determinar numero optimo de clusters

cols_num_norm_w <- na.omit(cols_num_norm)
fviz_nbclust(cols_num_norm_w, kmeans, method = "gap")

## Warning: did not converge in 10 iterations

## Warning: did not converge in 10 iterations

El numero optimo de clusters es 4

number_clusters <- 4

K-Means para clustering

km <- kmeans(cols_num_norm_w, centers = number_clusters, iter.max = 100)
km$size

## [1] 167 375 332 327

fviz_cluster(km, cols_num_norm_w)

Como se puede ver los datos si estan agrupados, él único que muestra un mínimo de dispersión es el primer cluster (de color rojo)

Clustering Jerarquico

hc<-hclust(data_dist, method = "ward.D2")
plot(hc, cex=0.5, axes=FALSE)
rect.hclust(hc,k=number_clusters)

También se puede ver que los clusters estan agrupados. En efecto VAT estaba en lo correcto.

4. Dividiendo los datos en entrenamiento y prueba

Dividimos el test de datos en un 75% para entrenamiento y un 25% para pruebas. Utilizamos el metodo de bootstraping porque nos asegura que la distribución de la data toma en cuenta la población total. Es decir, no tendremos un sesgo genrado por cómo la información está organizada.

set.seed(5)
expectedResult <- cleanData$salePrice
partition <- createDataPartition(y=expectedResult,
                                 p=.75,
                                 list=F)

trainingSet <- cleanData[partition,]
testingSet <- cleanData[-partition,]

5. Ingeniería de los Datos

Visualizamos la matriz de correlación para ver tendencias en la data.

correlations <- cor(cleanData[,c("neighborhood", "yearBuilt", "overallCondition", "overallQuality", "houseZone", "houseUtilities")], use="pairwise.complete.obs")
corrplot(correlations, method="circle", type="lower",  sig.level = 0.01, insig = "blank")

De esta figura obtenemos las siguientes conclusiones: - El año en el que se hizo la casa está correlacionado con el vecindario en el que se encuentra. - La condición de la casa tienen una correlación ligeramente negativa. Puede ser que mientras el año es menor, la condición actual reduce. - La zona donde se encuentra la casa también está ligeramente relacionada a la calidad de los materiales.

Continuamos con un segundo diagrama

correlations <- cor(cleanData[,c("neighborhood", "yearBuilt", "overallCondition", "overallQuality", "pool", "fireplace", "foundationMaterial", "exteriorCondition", "additionalFeature", "livingArea" )], use="pairwise.complete.obs")
corrplot(correlations, method="circle", type="lower",  sig.level = 0.01, insig = "blank")

Este diagrama nos da otra información de las variables: - Mientras mejor sea la condición del exterior, mejor es la condición en general de la casa. - la pisicna está ligeramente correlacionada con tener una fogata y con la condición del exterior de la casa. - El material usado para la base de la casa está relacionado con el año de construcción de esta. - El area de vivienda de la casa incrementa conforme la calidad de la construcción lo hace. - Una área de vivienda mayor permite tener más comodidades, como fogata o una piscina.

Si hacemos una gráfica de puntos tomando en cuenta algunas variables importantes podemos obtener más información:

pairs(~yearBuilt+overallQuality+overallCondition+livingArea,data=cleanData,
   main="Matriz de correlación")

- la calidad de las casas ha mejorado con el paso de los años. - El área de vivienda aumenta ligeramente con el paso del tiempo.

Tomando este análisis en cuenta, se utilizarán las siguientes variables como candidatas a mejores predicciones: - yearBuilt - el año de construcción tiene una alta influencia en la calidad de la casa, vecindario donde se encuentra, condición actual y área de vivienda. - livingArea - el area de vivienda es un factor imporante porque indica cuánto espacio puede ser aprovechado en la casa. No es la misma área que necesita una familia de 2 personas a una de 5. Además, un área mayor permite tener comodidades como piscina o fogata. - overallQuality - La calidad de la casa también es importante para el precio. Esto nos indica que el material usado es resistente. - overallCondition - La condición actual de la casa da una mejor impresión al mostrarla a posibles clientes. - pool - Una piscina aumenta el precio dado que el cliente debe tener presupuesto para mantenerla. - fireplace - en Iowa (lugar del dataset) hay inviernos con nieve, los clientes agradecerán esta comodidad

6. Asegurando resultados consistentes

set.seed(5)

7. Modelo univariable

Estimación de precio usando el area de vivienda

livingArea se toma como la variable independiente y se utiliza para hacer el modelo univariado de regresión lineal para predecir el precio de las casas.

singleVariableModel2 <- lm(salePrice~livingArea, data = trainingSet)

Caracteristicas del modelo univariado entrenado:

summary(singleVariableModel2)

## 
## Call:
## lm(formula = salePrice ~ livingArea, data = trainingSet)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -337959  -28898    -713   22328  243922 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 16969.633   5029.813   3.374 0.000767 ***
## livingArea    108.156      3.152  34.311  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 53750 on 1095 degrees of freedom
## Multiple R-squared:  0.5181, Adjusted R-squared:  0.5177 
## F-statistic:  1177 on 1 and 1095 DF,  p-value: < 2.2e-16

Los coeficientes en este modelo son significativos. Sin embargo \(R^2\) no es considerado alto, por tanto el modelo no logra explicar la variación en salePrice en su mayoria.

Ecuación de regresión: \(salePrice = 108.16livingArea + 1.696963\times 10^{4}\)

ggplot(data = trainingSet, mapping = aes(x = livingArea, y = salePrice)) +
geom_point(color = "green", size = 2) +
geom_smooth(method = "lm", se = TRUE) +
labs(title = "Precio de venta ~ Área de vivienda", x = "Área de vivienda", y = "Precio de venta") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5))

## `geom_smooth()` using formula = 'y ~ x'

Analisis de residuales

Ecuación para predecir el precio de venta para el conjunto de prueba.

predSP2<-predict(singleVariableModel2, newdata = testingSet)

Predicción

head(predSP2)

##        1        4        5        6        8       15 
## 201917.1 202674.1 254697.4 164278.6 243016.5 152489.6

length(predSP2)

## [1] 363

Valores de los residuos del modelo

head(singleVariableModel2$residuals)

##            2            3            7            9           10           11 
##  28037.00676  13363.06021 106813.44777 -78939.06315 -15454.06173     47.72457

Gráficos para analizar residuales

plot(singleVariableModel2)

En el gráfico Residuals vs Fitted se puede ver que los datos se encuentran alrededor de 0 pero no de forma aleatoria especificamente ya que la mayor cantidad de puntos se centran en un lugar en específico.

Chequeo de distribución aleatoria de los puntos:

hist(singleVariableModel2$residuals)

boxplot(singleVariableModel2$residuals)

qqnorm(singleVariableModel2$residuals)
qqline(singleVariableModel2$residuals, col="red")

Según los gráficos se puede observar que en el histograma obtenido no tiene una forma normal, ya que los datos se encuentran distribuidos hacia la derecha. Por otra parte en el gráfico q-q los extremos se alejan de la línea y en la caja de bigotes se observa que hay varios datos que se encuentran en los extremos.

library(nortest)
lillie.test(singleVariableModel2$residuals)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  singleVariableModel2$residuals
## D = 0.099646, p-value < 2.2e-16

Haciendo una prueba de normalidad el valor de p es menor que 0.05 por tanto la hipótesis nula de normalidad se rechaza y se afirma que los datos de los residuos no siguen una distribución normal.

Precio de venta para el conjunto de prueba

predSPP2<-predict(singleVariableModel2, newdata = testingSet[,c(19,21)])
library(caret)
RMSE(predSPP2,testingSet$salePrice)

## [1] 62571.91

plot(testingSet$salePrice,col="blue")
points(predSPP2, col="red")

summary(testingSet$salePrice-predSPP2)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -467188.0  -30742.2   -1108.3     178.9   21154.0  339005.7

En base a los resultados anteriores no se puede afirmar que se tenga una buena distribución ya que el RMSE es bastante alto, los datos en el gráfico se encuentran principalmente en la parte inferior del gráfico y no en diagonal. Además, en el resumen de estadisticas no se encuentran alrededor de 0 y tienen valores muy altos.

8. Modelo multivariable

Utilizamos todas las variables como una primera aproximación.

allVariablesModel <- lm(salePrice ~ ., data=trainingSet)
summary(allVariablesModel)

## 
## Call:
## lm(formula = salePrice ~ ., data = trainingSet)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -244657  -18333   -2015   15695  203644 
## 
## Coefficients: (2 not defined because of singularities)
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -5.440e+05  1.261e+05  -4.315 1.75e-05 ***
## neighborhood       -2.628e+03  2.782e+02  -9.449  < 2e-16 ***
## houseStyle         -2.144e+03  1.088e+03  -1.971 0.048994 *  
## houseZone           7.459e+03  2.336e+03   3.193 0.001449 ** 
## houseUtilities             NA         NA      NA       NA    
## roadAccess         -1.170e+03  1.185e+03  -0.987 0.323691    
## remodelated         7.964e+03  2.297e+03   3.467 0.000548 ***
## roofStyle           6.302e+03  1.252e+03   5.034 5.65e-07 ***
## roofMaterial        2.306e+02  1.878e+03   0.123 0.902286    
## overallQuality      1.594e+04  1.275e+03  12.501  < 2e-16 ***
## overallCondition    4.513e+03  1.112e+03   4.060 5.26e-05 ***
## exteriorCondition  -1.826e+03  3.056e+03  -0.597 0.550408    
## foundationMaterial  1.462e+03  2.112e+03   0.692 0.489044    
## basement            1.104e+04  6.783e+03   1.628 0.103876    
## basementCondition          NA         NA      NA       NA    
## fireplace           6.022e+03  2.403e+03   2.506 0.012351 *  
## garageArea          5.114e+01  6.185e+00   8.269 4.06e-16 ***
## pool                1.218e+05  1.916e+04   6.355 3.10e-10 ***
## additionalFeature   1.013e+03  5.901e+03   0.172 0.863726    
## livingArea          5.190e+01  2.877e+00  18.037  < 2e-16 ***
## yearBuilt           2.454e+02  6.453e+01   3.802 0.000152 ***
## saleType            4.877e+02  6.355e+02   0.767 0.442992    
## mSZoning           -9.137e+03  3.183e+03  -2.870 0.004182 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32540 on 1046 degrees of freedom
##   (30 observations deleted due to missingness)
## Multiple R-squared:  0.8264, Adjusted R-squared:  0.8231 
## F-statistic:   249 on 20 and 1046 DF,  p-value: < 2.2e-16

El valor de R cuadrado con este modelo es de 0.8231. Un valor arriba de 0.7 es considerado bueno para problemas de correlación. (Fernando, 2021) https://www.investopedia.com/terms/r/r-squared.asp#:~:text=In%20finance%2C%20an%20R%2DSquared,depend%20on%20the%20specific%20analysis.

Prediciendo valores con este modelo

predictionAllVariables <- predict(allVariablesModel, testingSet, type="response")

## Warning in predict.lm(allVariablesModel, testingSet, type = "response"):
## prediction from a rank-deficient fit may be misleading

allVariablesOutput <- cbind(testingSet, predictionAllVariables)

Observando primeros valores de la predicción

head(allVariablesOutput[,c("salePrice", "predictionAllVariables")])

##    salePrice predictionAllVariables
## 1     208500               216973.9
## 4     140000               214242.5
## 5     250000               299280.3
## 6     143000               159179.9
## 8     200000               224751.9
## 15    157000               160281.3

plot(allVariablesModel)

La gráfica de Residuales vs Fitted no presenta los valores distribuídos aleatoriamente de manera horizontal centrados en y=0. Sin embargo, tampoco está del todo mal si consideramos el valor de R cuadrado. La gráfica Q-Q también nos da información útil: paraece que los valores se aproximan a una línea recta entre los cuartiles teóricos -2 y 2.

complete_predicted <- na.omit(allVariablesOutput)
rmseAllVariables <- rmse(complete_predicted$salePrice,complete_predicted$prediction)

El valor del RMSE es 49481.26.

9. Analizando el modelo

Graficamos la matríz de correlación para todas las variables del modelo.

modelCorrelations <- cor(allVariablesOutput, method="pearson")
modelCorrelations[is.na(modelCorrelations)] <- 0

corrplot(modelCorrelations, type="upper", order="hclust", tl.col="black", tl.srt=45)

En el resumen generado arriba se puede observar que las variables más significativas son: neighborhood, remodelated, roofStyle, overallQuality, overallCondition, garageArea, pool, livingArea, yearBuilt.

La gráfica de los valores ajustados vs los residuales nos indica que el modelo puede tener un problema de overfitting. Es decir, no logrará generalizar correctamente los valores del set te prueba.

10. Modelo con variables significativas

newMultivariableModel <- lm(salePrice ~ neighborhood + remodelated + roofStyle + overallQuality + overallCondition + garageArea + livingArea + yearBuilt, data=trainingSet)
summary(newMultivariableModel)

## 
## Call:
## lm(formula = salePrice ~ neighborhood + remodelated + roofStyle + 
##     overallQuality + overallCondition + garageArea + livingArea + 
##     yearBuilt, data = trainingSet)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -250231  -18907   -2185   15338  281048 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -6.372e+05  1.135e+05  -5.612 2.53e-08 ***
## neighborhood     -2.728e+03  2.539e+02 -10.743  < 2e-16 ***
## remodelated       7.881e+03  2.311e+03   3.410 0.000673 ***
## roofStyle         7.133e+03  1.243e+03   5.740 1.23e-08 ***
## overallQuality    1.607e+04  1.223e+03  13.145  < 2e-16 ***
## overallCondition  5.112e+03  1.011e+03   5.055 5.06e-07 ***
## garageArea        5.169e+01  6.167e+00   8.382  < 2e-16 ***
## livingArea        5.341e+01  2.617e+00  20.403  < 2e-16 ***
## yearBuilt         3.046e+02  5.663e+01   5.379 9.19e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33470 on 1088 degrees of freedom
## Multiple R-squared:  0.8143, Adjusted R-squared:  0.813 
## F-statistic: 596.4 on 8 and 1088 DF,  p-value: < 2.2e-16

plot(newMultivariableModel)

La reducción de variables disminuyó el sobreajuste. Sin embargo, también disminuyó el valor de R cuadrado a 0.81. Hay que tomar en cuenta que esto es para los datos de entrenamiento, así que pondremos a prueba al modelo:

newModelPrediction <- predict(newMultivariableModel, testingSet, type="response")
newModelOutput <- cbind(testingSet, newModelPrediction)

Observando primeros valores de la predicción

head(newModelOutput[,c("salePrice", "newModelPrediction")])

##    salePrice newModelPrediction
## 1     208500           220360.8
## 4     140000           209398.8
## 5     250000           298290.8
## 6     143000           157311.9
## 8     200000           225137.8
## 15    157000           154551.2

completeNewModelPredicted <- na.omit(newModelOutput)
rmseNewModel <- rmse(completeNewModelPredicted$salePrice,completeNewModelPredicted$newModelPrediction)

El valor del RMSE es 45207.99. Lo cuál es menor al valor del modelo anterior (49481.26). Esto es un indicador que el nuevo modelo es mejor generalizando la información.

11. Modelos con conjunto de pruebas

Modelo 1 - Univariable

svmodel7 <- lm(salePrice~livingArea, data = testingSet)

Caracteristicas del modelo univariado entrenado:

summary(svmodel7)

## 
## Call:
## lm(formula = salePrice ~ livingArea, data = testingSet)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -452159  -31302   -2239   20302  341908 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 22812.106   9650.226   2.364   0.0186 *  
## livingArea    104.457      5.925  17.629   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 62710 on 361 degrees of freedom
## Multiple R-squared:  0.4626, Adjusted R-squared:  0.4611 
## F-statistic: 310.8 on 1 and 361 DF,  p-value: < 2.2e-16

Podemos ver que \(R^2\) (0.46) es más bajo que cuando se utilizó el conjunto de entrenamiento (\(R^2\)) por tanto la eficiencia del algoritmo para predecir el precio de las casas es menor y la regresión lineal se acopla de menor forma al conjunto de datos.

Modelo 2 - Multivariable

allVM8 <- lm(salePrice ~ ., data=testingSet)
summary(allVM8)

## 
## Call:
## lm(formula = salePrice ~ ., data = testingSet)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -227319  -20567   -3701   16035  256371 
## 
## Coefficients: (1 not defined because of singularities)
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -1.089e+06  3.052e+05  -3.570 0.000411 ***
## neighborhood       -2.951e+03  5.956e+02  -4.955 1.16e-06 ***
## houseStyle         -3.980e+03  2.489e+03  -1.599 0.110793    
## houseZone           1.043e+04  5.596e+03   1.864 0.063216 .  
## houseUtilities     -1.325e+04  2.349e+04  -0.564 0.573130    
## roadAccess          2.003e+03  2.855e+03   0.702 0.483428    
## remodelated         1.531e+04  5.531e+03   2.768 0.005962 ** 
## roofStyle           7.473e+03  2.927e+03   2.553 0.011134 *  
## roofMaterial        1.381e+04  3.430e+03   4.025 7.09e-05 ***
## overallQuality      1.611e+04  3.139e+03   5.130 4.95e-07 ***
## overallCondition    2.801e+03  2.778e+03   1.009 0.313937    
## exteriorCondition   5.398e+03  8.289e+03   0.651 0.515399    
## foundationMaterial -2.176e+03  4.350e+03  -0.500 0.617261    
## basement            7.320e+03  1.896e+04   0.386 0.699684    
## basementCondition          NA         NA      NA       NA    
## fireplace           5.562e+03  5.428e+03   1.025 0.306254    
## garageArea          2.862e+01  1.542e+01   1.856 0.064295 .  
## pool               -8.657e+04  2.407e+04  -3.596 0.000372 ***
## additionalFeature   4.398e+02  1.038e+04   0.042 0.966226    
## livingArea          5.401e+01  6.367e+00   8.483 7.53e-16 ***
## yearBuilt           5.036e+02  1.547e+02   3.255 0.001251 ** 
## saleType           -1.773e+03  1.933e+03  -0.917 0.359748    
## mSZoning           -6.299e+03  7.541e+03  -0.835 0.404122    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 43040 on 329 degrees of freedom
##   (12 observations deleted due to missingness)
## Multiple R-squared:  0.7533, Adjusted R-squared:  0.7376 
## F-statistic: 47.85 on 21 and 329 DF,  p-value: < 2.2e-16

El valor de \(R^2\) con este modelo es de 0.75 y aunque es menor a cuando se utilzó el conjunto de entrenamiento (\(R^2\) = 0.8231) Aún es considerado aceptable para problemas de correlación. (Fernando, 2021) https://www.investopedia.com/terms/r/r-squared.asp#:~:text=In%20finance%2C%20an%20R%2DSquared,depend%20on%20the%20specific%20analysis.

Modelo 3 - Multivariable con variables significativas

allVM10 <- lm(salePrice ~ neighborhood + remodelated + roofStyle + overallQuality + overallCondition + garageArea + livingArea + yearBuilt, data=testingSet)
summary(allVM10)

## 
## Call:
## lm(formula = salePrice ~ neighborhood + remodelated + roofStyle + 
##     overallQuality + overallCondition + garageArea + livingArea + 
##     yearBuilt, data = testingSet)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -323417  -21063   -4141   17198  293648 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      -7.713e+05  2.589e+05  -2.980  0.00308 ** 
## neighborhood     -3.080e+03  5.396e+02  -5.708 2.42e-08 ***
## remodelated       1.581e+04  5.598e+03   2.825  0.00500 ** 
## roofStyle         7.315e+03  2.912e+03   2.512  0.01245 *  
## overallQuality    1.908e+04  2.857e+03   6.678 9.38e-11 ***
## overallCondition  5.168e+03  2.434e+03   2.124  0.03438 *  
## garageArea        3.123e+01  1.462e+01   2.136  0.03334 *  
## livingArea        5.017e+01  5.928e+00   8.463 7.01e-16 ***
## yearBuilt         3.716e+02  1.285e+02   2.892  0.00407 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 45310 on 354 degrees of freedom
## Multiple R-squared:  0.7249, Adjusted R-squared:  0.7186 
## F-statistic: 116.6 on 8 and 354 DF,  p-value: < 2.2e-16

En este caso a pesar de \(R^2\) (0.72) aún es considerado aceptable para problemas de correlación (Fernando, 2021) su valor fue reducido considerablemente vs si se hubiera hecho con el conjunto de datos de entrenamiento (0.81). https://www.investopedia.com/terms/r/r-squared.asp#:~:text=In%20finance%2C%20an%20R%2DSquared,depend%20on%20the%20specific%20analysis.

Es posible concluir que para los modelos tanto univariables como multivariables la eficiencia del modelo se ve reducida y un menor \(R^2\) nos indica que el modelo que el modelo se ajusta menos a los datos y por tanto hay una menor precisión en la predicción de los datos.

12. Efectividad de modelos

modelo	rCuadrado
Lineal	0.5200
Multivariable	0.8231
Multivariable Significativas	0.8130

El valor de R cuadrado pareciera ser mejor en el modelo Multivariable que usa todas las variables numérics. Sin embargo, hay que tomar en cuenta que esta medida es para el set de entrenamiento. Si se observa el análisis de residuales para ese modelo se ve que este está sobreajustado.

Graficando Modelo Univariable (con LivingArea)

Graficando Modelo Multivariable

En azul tenemos el precio correcto y en rojo los valores predecidos para ese mismo conjunto de datos. Podemos observar que el modelo no es muy bueno generalizando precios muy altos. Lo cuál confirma la forma que se ve en la gráfica de residuales.

Graficando Modelo Multivariable con variables significativas

Este modelo a pesar de ser parecido al anterior puede generaizar mejor los valores de prueba. Esto se ve claramente en las observaciones con precios bajos.

Por otro lado, al comparar el RMSE de los modelos (62571.91, 49481.26 y 45207.99), obtenemos que el tercero tiene un menor valor, es decir, es un indicador que es mejor generalizando los resultados.

En conclusión, el modelo que mejor generaliza los precios de las casas es el modelo multivariable que utiliza únicamente las variables importantes: neighborhood, remodelated, roofStyle, overallQuality, overallCondition, garageArea, pool, livingArea, yearBuilt. Hay que tomar en cuenta que podríamos tener mejores resultados utilizando un modelo polinómico en vez de lineal porque algunas características, como livingArea no crecen de forma lineal.

notebook

2023-03-02

Hoja de trabajo 03 - Regresión Lineal

2. Análisis Exploratorio

¿Cuál es la Zona más cara?

¿En qué vecindario se ubican las casas con mayor área en su terreno?

¿En qué año hubo más remodelaciones?

¿Cuáles son los materiales más utilizados para la fundación de las casas?

¿Qué tipos de casa son más comunes?

¿Que vecindarios tienen un mayor rango de precios?

¿Cuál proporción de casas en cada zona?

¿Cuáles son las características más comunes de las casas que se venden por encima del precio medio?

¿Cuántas casas tienen piscina o garaje? ¿Como se relacionan estas características con el precio de venta?

¿Cómo se distribuyen las ventas de casas a lo largo del tiempo entre los diferentes periodos?

Preprocesamiento de datos

Borrando valores inneceesarios

3. Clustering

Hopkins

Evaluación Visual de Tendencia de la data (VAT)

Determinar numero optimo de clusters

K-Means para clustering

Clustering Jerarquico

4. Dividiendo los datos en entrenamiento y prueba

5. Ingeniería de los Datos

6. Asegurando resultados consistentes

7. Modelo univariable

Estimación de precio usando el area de vivienda

Caracteristicas del modelo univariado entrenado:

Analisis de residuales

Predicción

Valores de los residuos del modelo

Gráficos para analizar residuales

Chequeo de distribución aleatoria de los puntos:

Precio de venta para el conjunto de prueba

8. Modelo multivariable

9. Analizando el modelo

10. Modelo con variables significativas

11. Modelos con conjunto de pruebas

Modelo 1 - Univariable

Caracteristicas del modelo univariado entrenado:

Modelo 2 - Multivariable

Modelo 3 - Multivariable con variables significativas

12. Efectividad de modelos

Graficando Modelo Univariable (con LivingArea)

Graficando Modelo Multivariable

Graficando Modelo Multivariable con variables significativas